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We study the effects of finite size and of vacancies on tlie pliotonic band gap recently predicted 
for an atomic diamond lattice. Close to a Jg = —>■ Je = 1 atomic transition, and for atomic 
lattices containing up to Af « 3 x 10^ atoms, we show how the density of states can be affected 
by both the shape of the system and the possible presence of a fraction of unoccupied lattice sites. 
We numerically predict and theoretically explain the presence of shape-induced border states and 
of vacancy- induced localized states appearing in the gap. We also investigate the penetration depth 
of the electromagnetic field which we compare to the case of an infinite system. 
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I. INTRODUCTION 



That of waves propagation in periodic potentials con- 
stitute a problem shared by several domains of classical 
and quantum physics, ranging from the study of electron 
motion in metals [ij, to that of X- and 7-ray scatter- 



ing by crystals [iHSL, and of light by photonics crystals 
and metamaterials ^. Periodicity leads to an organiza- 
tion of modes according to bands, and to the possible 
presence of band gaps, i.e. energy intervals where modes 
are absent. A periodic system is by definition infinitely 
extended, hence not physical. Nonetheless, predictions 
made on the base of infinite systems can become really 
satisfactory for systems large enough, as in solid-state 
physics, and present the advantage to benefit from the 
Bloch theorem, and to be solved in the reciprocal space 
avoiding typical real space oscillating functions. Mod- 
els based on infinite systems may however present some 
subtleties related to the way in which the infinite limit- 
ing process is performed, often requiring ad hoc Ewald's 
summations type strategies. 

The recent experimental realization of a Mott phase 
with ultracold atomic gases [a, Qi i-^- of artificial crys- 
tals made by single atoms trapped at the nodes of laser 
optical lattices, leads to the necessity of understanding 
the features of the band structure of light interacting 
with such systems. The peculiarity of this new system 
is that it presents several remarkable features: incident 
light scatters on point-like elementary quantum objects 
with an internal energy level structure, and a quantum 
delocalized position in space [7[; the lattice periodicity 
is of the order of the incident light wavelength, allowing 
the exploration of the entire Brillouin zone and hence of 
possible band gaps |8| ; experiments reached a remarkable 
accuracy and control, permitting the realization of ultra- 
precise atomic clocks j9l-[ll|. First attempts toward the 
description of such a system overlooked divergence prob- 
lems, resulting in non correct prediction of band gaps 
p^ Il3j , or were based on a ad hoc ultraviolet regularisa- 



tion procedure allowing to explore only a particular class 
of lattice geometries not presenting any band gap [14| . 
Photonic band gaps of ID cold atomic vapors have been 
realized [l5|, and exploited to generate o ptic al paramet- 
ric oscillation with distributed feedback [l^. Scattered 
photons have been suggested as a signature of the Mott 
insulator and superfluid quantum states [l7| , and studied 
in the framework of polaritons [l8| , excitons and cavity 
polaritons [l9| . Recently, by exploiting a microscopic 
theory of light-atom interaction [20[, and by explicitly 
introducing the presence of the unavoidable atomic quan- 
tum motion, it was possible to naturally regularize the 
divergences in a way independent of the lattice geometry, 
and at the same time to study the quantitative effects of 
the quantum atomic motion on the band structure [7|. 
The explicit dependence of the photonic band structure 
on quantum features, as the atomic internal energy lev- 
els and the external atomic quantum motion, allows to 
consider this artificial structure as an example of quan- 
tum metamaterial |2l| . In the framework of the Fano- 
Hopfield self-consistent quadratic theory [3, [l^ [g^], it 
was also possible to find an exact solution valid for the 
full Brillouin zone and for arbitrary Bravais and non- 
Bravais lattices, allowing the prediction of the diamond 
as the first 3D atomic lattice geometry presenting a com- 
plete photonic band gap Q [S^]- Further investigations 
suggested to add external magnetic fields to open band 
gaps in other geometric structures [2a |. 

In cold atom realizations of 3D optical lattices, the 
atomic Mott state extends over 10 — 20 lattice sites, so 
a natural question regards the features of the band gap 
in this finite size system. A further question concerns 
the effects of an imperfect finite portion of a lattice con- 
taining site defects, i.e. a fraction of vacancies resulting 
in a not complete filling of the lattice. The experimental 
interest of these issues is related to the fact that both the 
finite size and vacancy effects, separately, could in prin- 
ciple drastically affect the presence and the experimental 
visibility of the band gap due to the appearance of states 
in the gap region. The main questions we address in 



this paper arc: What docs happen to the band gap for 
systems of reahstic sizes and of different shapes? What 
is the fraction of vacancies which still permit to have 
a reasonable band gap visibility? What is the value of 
the penetration depth of an electromagnetic wave in the 
atomic diamond lattice for finite and infinite systems, i.e. 
how is it affected by finite size effects? Even if we discuss 
in detail the case of a diamond lattice, we will present a 
general formulation and will discuss main features which 
will remain valid for other lattice geometries. 

The paper is organized as follows. In section|ll]we illus- 
trate the model we use, and the resulting main equations 
we solve. In section Hill we present and discuss a numeri- 
cal study on the density of states and on the penetration 
length in a finite size system, possibly in presence of im- 
perfections due to vacant sites in the lattice. In section 
II VI we provide an analytical analysis to support and illus- 
trate some of the main features of the numerical findings. 
We conclude in section IVl 



II. THE MODEL 

We consider a system made by a collection of N iden- 
tical atoms having fixed positions and an optical dipolar 
transition between a Jg = electronic ground state and 
a Je = 1 electronic excited state ^26| . Such a transition is 
available in appropriate atomic species, such as strontium 
where it was already used to study coherent propagation 
of light in an atomic ensemble [24|- In our model, the 
atomic dipoles are coupled by the electromagnetic field 
they radiate, and in the regime of low atomic excitations, 
the resulting eigenmodes of the mean atomic dipoles are 
given by the solutions of the eigenvalue problem [ij, [20] 



j=l p=x,y,z 

h{uj - ij)di^a- (1) 



Here di_a is the component along the direction a = x,y, z 
of the mean dipole carried by the atom i, u — ij is the 
mode eigenfrcquency (it is complex in general with 7 > 0, 
and may be measured as suggested in Q), wq and F 
are the single atom resonance frequency and spontaneous 
emission rate. The tensor gaj3{r) gives the a component 
of the electric field at the position r radiated by a dipole 
oscillating along the direction /? at the origin of coordi- 
nates, Ea{r) = —gapi'r)dj3/d'^, d being the atomic dipole 
moment such that F = d'^uiQ / {STreohc^) . Here we con- 
sider the case where a; — 17 is very close to ujq-, so that 
gap can be evaluated for a dipole oscillating at the reso- 
nance frequency; introducing the vacuum wavenumber 



we take 
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Our first expression in ([3]) for ,90/3 (r) differs by a scalar 
6(r) contribution from the usual expression for the elec- 
tric field radiated by a dipole, see Eqs. (4.20,9.18) of |28| : 
this ensures compatibility with our previous works and it 
is of course irrelevant here since atoms are never at the 
same position |33|. The first expression is particularly 
useful to directly extract its Fourier transform, needed 
in the Bloch-description of infinite systems (see section 
IIVP , while the second one (which differs from the first one 
by another scalar S{r) contribution) has the well know 
dipole-dipole interaction form, and will be used in nu- 
merical calculations on finite-size systems in section Hill 
Equation ([T]) allows one to determine the density of 
states of the system. In case an infinite number of atoms 
are periodically arranged at the nodes of a diamond lat- 
tice, it has been shown that the system may exhibit an 
omnidirectional photonic band gap [8| . Here, by numer- 
ical solution of ([H we investigate the fate of such a gap, 
in situations close to realistic experimental ones, where 
the number of atoms is finite and/or there are vacan- 
cies in the lattice. A further interesting quantity related 
to the occurrence of a gap is the so-called "penetration 
depth" ^: an incident electromagnetic wave at a fre- 
quency in the band gap cannot penetrate the medium, 
and its amplitude will decay exponentially over a char- 
acteristic distance ^. In order to calculate such a length 
we consider a point-like dipolar source immersed in the 
atomic medium, and we extract ^ from the total field and 
the induced dipole spatial profiles: we fix at the position 
Vs a forced dipole d^ = d^ e~*"=*, the atomic dipoles at 
the positions r^ will reach a steady state di^a = rfi.a e~*"=* 
given by the linear system 



hiiJs ^ ^0) + i 



hT 



N 



j = l P=x,y,z 

- J2 3«/3(l■^-r.)d^. (4) 



III. 



=x,y,z 



NUMERICAL RESULTS FOR A FINITE 
SIZE SYSTEM 



Wo 

c 



(2) 



In this section we study a system of N atoms at the 
nodes of a diamond lattice. We recall that the diamond 
lattice is formed by the superposition of two copies of the 



same Bravais lattice: the fee lattice of lattice constant a, 
generated by the three basis vectors 

ei = (0, a/2, a/2), 62 = (a/2,0, a/2), eg = (a/2, a/2,0), 

(5) 
and a second fee lattice obtained by translating the first 
lattice by the vector (a/4, a/4, a/4). The corresponding 
basis of the reciprocal lattice is 

§1 = (— 27r/a, 27r/a, 27r/a), 62 = (27r/a, — 27r/a, 27r/a), 
§3 = (27r/a, 27r/a, -27r/a). (6) 

In our simulations, the atoms occupy a finite region in 
space, which can be a ball or a cube centered at the 
origin of the coordinates. From the numerical solution 
of ([T]) we extract the density of states for the case of a 
unit filling factor (section IIII A[) and for the case with 
a low concentration of vacancies (section |IIIB[) . Finally, 
we analyze the penetration depth in section llll CI solving 



A. Finite size effects on the density of states 

In this section we discuss the density of states obtained 
by solving equation ([1]) for a finite size diamond lattice, 
in the absence of vacancies. In particular, in Fig. [T] we 
show the density of states p{u!) for a number of atoms 
corresponding to typical experimental realizations N « 
2.5 X 10^. Here p{uj) is defined as the distribution of 
the real part w of the complex spectrum of equation ([l} , 
normalized as J p{uj)dLj = 6/Vl, where Vl = a'^/4 is the 
volume of the direct lattice primitive cell, in order to 
facilitate the comparison with the infinite system results 
of 0, plotted as a bar histogram in the figure. If the 
atoms occupy a ball (see the black solid line) we observe 
partial filling of the spectral gap, most pronounced in 
the upper region. On the contrary the region close to the 
lower border of the gap remains relatively weakly affected 
by the finite size of the system, considering the sharp rise 
of p{uj) to the left of this border. The remaining part of 
the density of states remains very close to the one of 
the infinite system. If the atoms occupy a cube (see the 
red solid line) the finite size effects are quite different. 
Two peaks appear, a very pronounced one in the middle 
of the gap (at (w — wo)/r w —3.2), and a second one 
(at (oj — a;o)/r w 0.5). We investigated the nature of 
the states belonging to the peak in the gap, by looking 
at 10 successive eigenstates of ([1]), finding that they are 
"border states" : they reach their maxima in a spherical 
shell of radius « 5a, and decay exponentially towards the 
center of the cube with a law 



|d,p= Y. \d,^a\'< 



-22+4.6ri/a 



(7) 




where the dipole eigenvectors are normalized to the maxi- 
mum value of their modulus equal to unity. This suggests 



{(0-%W 



Figure 1: (Color online) Density of the real part of the eigen- 
frequencies p{uj) obtained from Eq.([l}, in the absence of va- 
cancies and for koa = 2 where a is the fee lattice constant. 
Red solid lines: finite system with a cubic shape of side of 
length 14a, and A'' — 2.7 x 10'*. Black solid lines: finite system 
with a spherical shape of diameter 18a, and A^ = 2.4 x 10''. 
The histogram provides the same quantity for an infinite sys- 
tem [3. Each of the three curves is composed of 250 bins. Vl 
is the volume of the direct lattice primitive cell. The inset is 
a magnification. 



a value of the penetration depth of the order of 0.5a, in 
agreement with the calculation done in section llll CI 

In Fig. [5] we show the distribution of the eigenvalues of 
Eq.(IT]) in the complex plane, restricted to small values of 
7/r. In this region, the figure shows that the band gap 
is not filled, apart from two narrow intervals of values of 
w, in the center and close to the upper border of the gap. 
Then, in the finite size system, the partial filling of the 
gap is mostly due to eigenvalues with larger values of 7/r. 
The smallest values of 7/r we obtained are « 2 x 10~^. 
The real part of the corresponding eigenvalues are located 
on the borders of the band gap for the infinite system, 
marked in the figure by vertical dashed lines, and on the 
upper bound of the values of ui represented in Fig. [2] for 
the finite system, i.e. around {uj — uj{))/T ~ 9.5. 



B. Effects of vacancies on the density of states 

In this section we address the case where the finite 
size diamond lattice in not perfectly filled, presenting a 
concentration 1 — p of defects made by the presence of a 
random uniform distribution of not-occupied lattice sites. 
In Fig.[3]we show the density of states /c(w) obtained solv- 
ing Eq. ([1]) for atoms occupying a ball, as a function of 
the lattice filling factor p. The figure, and its inset, show 
that already a small concentration of vacancies equal to 
l—p~ 0.99 (red solid line) produces a remarkable signa- 
ture in the density of states manifested by the appearance 
of a pronounced peak in the middle of the band gap, at 
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Figure 2: (Color online) Complex eigenvalues cj — iy obtained 
from Eq.(IT|. The system is of finite size, in the absence of 
vacancies, for koa = 2, with a spherical shape of diameter 
18a, and A'^ = 2.4 x IC* atoms. The two vertical red dashed 
lines give the borders of the band gap of the infinite periodic 
system. 



(w — CiJo)/r « —3.08. We explain the nature of the peak 
with the presence of single-vacancy states localized at 
the vacancy position. Since the vacancy concentration 
is small, most frequent vacancy states have a single-site 
nature. In section IIVCI we theoretically calculate the 
value of the single-vacancy state frequency, signaled in 
the inset by a black vertical dotted line, which seems 
to coincide quite satisfactorily with that of the numeri- 
cally observed peak. By increasing the vacancy concen- 
tration. Fig. [3]shows for 1 — p = 0.05 the occurrence of a 
clear second peak in the gap, which seems to match quite 
well the frequency of a two-vacancy in-gap state calcu- 
lated in section HV C[ see the red vertical dotted line at 
(a; — a;o)/r ~ —4. Peaks corresponding to other two- 
vacancy states predicted in section IIVCI are less visible 
(see the other vertical dotted lines in the inset of Fig. [3]) . 
Further increase of the concentration of vacancies pro- 
duces a gradual filling of the band gap, whose visibil- 
ity completely deteriorates for a vacancy concentration 
around 1 — p= 0.2. 

In Fig. m for exactly the same spherical system with 
a vacancy concentration of 1 — p = 0.2, we show the 
distribution of the eigenvalues of Eq.([l]) in the complex 
plane, restricted to small values of 7/r. The figure shows 
that the band gap is completely filled. The states filling 
the gap, for such a large vacancy concentration, are com- 
pletely delocalized over the entire system size, and have 
a spectral imaginary part mostly concentrated around 
7/r sa 10-2, with 7/r > 10-3. 

In Fig. Owe study the effect of vacancies on a system 
of cubic shape. The figure shows that for a concentration 
of vacancies 1 — p = 0.01 (red sohd line) two peaks are 
present in the band gap. They have a different origin: 
the first one, that at smallest energy, in nothing but the 
peak related to shape-induced states, already present in 




(a>-a)o)/r 



Figure 3: (Color online) Density of the real part of the eigen- 
frequencies p{io) obtained from Eq.fT]), for different concentra- 
tions of vacancies, that is for various filling factors p, and for 
koa — 2. The finite system has a spherical shape of diameter 
18a, and A'' = 2.4 x 10 for p = 1. The histogram provides the 
same quantity for an infinite system with no vacancies [g| . Vl 
is the volume of the direct lattice primitive cell. The inset is 
a magnification, where the vertical dotted lines correspond to 
frequencies of the single vacancy in-gap state (black, central) 
and to two- vacancy in-gap states {R2 — Ri = ei, jli — 1I2 — 1 
in red, outer; R2 — Ri = ae^, fli = 2, fi2 = 1 in blue, inner; 
these quantities are defined in appendix |B]) theoretically pre- 
dicted in section IIVCI Decreasing values of p correspond to 
increasing values of p{uj) in the band gap. 



the absence of vacancies (see black solid line, and the 
discussion in section IIII Ap . The second peak is instead 
the signature of single-vacancy localized states, and its 
position is the same of that shown in Fig. |3] for spherical 
shape at the same vacancy concentration. 



C. Penetration depth 

To numerically calculate the penetration depth ^ for 
a diamond finite-size atomic lattice we numerically solve 
the forced dipole equation ^ for a point-like dipolar os- 
cillating source d^ = d^e-'"'^* at the position Vg (ap- 
proximately at the center of the system), and with ujs in 
the band gap. Solutions of Eq.(|3|) provide the induced 
atomic dipoles amplitudes di^a at the lattice positions r^. 

We extract ^ according to different methods. The first 
method is based on the direct analysis of the induced 

dipoles, and consist in averaging the norm \/X)q Mi,Qp 
on spherical shells of radius « u = ||r — rs|| centered 
at the source position. We then obtain an average real 
dipole function ^{u) that we fit in a certain range of u 
(where the behavior of d{u) is clearly exponential over 
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Figure 4: (Color online) Complex eigenvalues cj — i^ obtained 
from Eq.([l}. The system is of finite size, in presence of va- 
cancies, that is with a filling factor p = 0.8, for koa — 2, with 
a spherical shape of diameter 18a, and N = 1.9 x 10* atoms. 
The two vertical red dashed lines give the borders of the band 
gap of the infinite periodic system. 



method provides the results presented by red squares in 
Fig. [5] Its specialisation to the analysis of the penetra- 
tion depth along some given direction (without averaging 
over spherical shells) is straightforward, and leads to the 
filled diamonds and circles in Fig. [Tji and b, respectively. 
The second method is based on the calculation of the 
total electric field amplitude generated by the source and 
induced dipoles obtained by (H)) : 



u ^ d 

Ea{r) = -y^5a/3(r-rs)--| -^^ga^(r-r,)-^, 



'd2 



(9) 
We evaluate Ea{r) on three lines, parallel to the Carte- 
sian axes and passing trough the source position Vg ■ We 



first average the norm \/J2a l-^a(r)p on the two direc- 
tions (±) of the three axes a, then we obtain and fit the 
six corresponding average real electric functions £a (u) 
as 



£t\u)^Kt' 



=-«/«i 



(10) 
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Figure 5: (Color online) Density of the real part of the eigen- 
frequencies p(a;) obtained from Eq.(IT]), for two concentrations 
of vacancies, that is for the filling factors p — 1 and p = 0.99, 
and for koa = 2. The finite system has a cubic shape of side 
14a, and A'^ — 2.7 x lO** for p = 1. The histogram provides the 
same quantity for an infinite system with no vacancies [g| . Vl 
is the volume of the direct lattice primitive cell. We note the 
double peak structure for p = 0.99 (see text). The vertical 
dotted line corresponds to the frequency of the single vacancy 
in-gap state theoretically predicted in section HV CI 



several decades) as 



V{u) = C 



-u/i 



(8) 



where ^ and C are the two fitting parameters. The fac- 
tor 1/m in dS]) is introduced to take into account the di- 
rect effect of the source which is dominant at small dis- 
tances, allowing to fit the function on a larger range. This 



obtaining six values of ^q, , whose average is presented 
by empty black circles in Figs. |6] and [Tja. 

In Fig. [6l it is apparent that the extractions of the 
penetration depth from Eq. ([5]) and from Eq. pU)) give 
different values. This shows that ^ is not isotropic, it 
depends on the considered direction of space, a prop- 
erty that will be recovered analytically in section IIVBI 
Whereas use of Eq. ([T0|) is expected to give the pen- 
etration depth along x axis, the first method, when it 
involves a directional average as in Eq. ([8]), is expected 
to pull out the maximal penetration depth (maximized 
over the directions of space). A second property, appar- 
ent in Fig. [6^, is the divergence of ^ at the borders of 
the infinite-medium forbidden gap (represented by verti- 
cal dashed lines at frequencies ui-mi, Wgup)- Fig. ^p even 
suggests that k vanishes there with a vertical slope. We 
indeed find that k^ vanishes linearly with ujs (not shown) , 
as also predicted analytically in section HVBI By a linear 
extrapolation of k^ as a function of w^, we get for the 
borders of the forbidden bands: 



Eq. H 



(-4.747, -1.948)(12) 



which are indeed quite close to the infinite medium re- 
sults 11: 



t^inf — "^0 l^i 



sup 



Wo 



r 



)~ (-4.743, -T962). 



(13) 



To better put in evidence the vanishing of k at the band 
edges, and to more easily compare the various methods. 



we show K as a function of (w^ — Winf)/(w, 



sup 



s) 



Fig. [3 with the band edges Wjnf and Wg,,,, deduced for 



^sup 



the finite-size simulations by linear extrapolation of 




(a^cOo)/r 



Figure 6: (Color online) Penetration depth ^ in (a) and its in- 
verse K in (b), as functions of the dipole source frequency ujs- 
Symbols (the lines are a guide to the eye) correspond to the 
numerical solution of Eq.(|4]) for a finite system of spherical 
shape, diameter 18a, filling factor p = 1, koa = 2, and con- 
taining N = 2.4 X 10* atoms. Red squares and black circles 
correspond to values obtained using the methods of Eq.® 
and of Eg. HlOp . respectively. The vertical dashed lines cor- 
responds to the borders (|13p of the band gap for the infinite 
periodic system. 
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Figure 7: Inverse penetration depth k = 1/^ along direction 
ei in (a) and direction e^; in (b), as a function of the source 
frequency u)s expressed through a change of variable mapping 
the band gap [a;inf,a;sup] onto R+. Same physical parameters 
as in Fig. [6l Filled diamonds in (a) and filled circles in (b): 
finite size system with first extraction method; empty circles 
in (b): finite size system with second extraction method; for 
those data, cjinf and Wsup were obtained by linear extrapo- 
lation of K^ to zero. Stars: for the infinite medium from a 



numerical evaluation of Eq. (|30p . Dashed lines: analytical 
predictions, close to the band borders, deduced from Eq. (|42p 
(see text). Note that the x and y axes are in logjQ and logj 
scale. 



IV. THEORY FOR THE INFINITE SYSTEM 



[3J|. This change of variable on Ug has the advantage 
of mapping the band edges to and 4-oo, respectively, 
which is then combined with a log-scale representation on 
both figure axes. This figure was produced for two par- 
ticular directions of penetration, along the direct lattice 
basis vector ei in Fig. [7]i, and along the Cartesian axis 
direction e^ in Fig. [7}d. First, in Fig. [TJj, it appears that 
the two extraction methods for the penetration depth in 
the finite-size system (the first method from the dipoles, 
see the filled circles; the second method from the electric 
field, see the empty circles) give compatible results if they 
are applied along the same direction (here e^;, which is 
equivalent to By or e^ due to symmetry of the diamond 
lattice). Second, in Fig. [7^ and b, the results of the finite- 
size systems are compatible with the ones (stars) for the 
infinite system in section IIVB[ and even if they do not 
cover a as large range for k, they do nicely follow the 
analytical prediction (dashed lines) for the vanishing of 
K close to the band edges. 



We show in this section that several features of the 
numerical simulations, such as the sharp rise of ^ close 
to the band gap borders and some peaks induced by va- 
cancies in p(w), can be interpreted analytically for an 
infinite system. In this case, a reformulation of (|ll4p 
in Fourier space is more appropriate. It is known how- 
ever that the resulting series over the reciprocal lattice 
present subtle convergence issues [Ijl that were over- 
looked in [l3, [13 ■ These issues were solved in Q by 
coupling each atomic dipole to a spatially smoothed ver- 
sion E^ (r) = / d^uE± (r — u)x(u) of the transverse elec- 
tromagnetic field operator Ej^(r), where the smoothing 
function x(u) may be taken as a positive rotationally in- 
variant function of unit integral and of small width b. 
This cuts off the dipolar coupling at high wavcnumber 
field modes and regularizes the theory for the infinite 
system. 

One then finds that two changes have to be applied to 
Eqs. (|1I4[1 . First, the function ^^^(r) has to be replaced 



by the smoothed function gapi^) such that 

gafsiri-Tj) ^ £u^ (fujga[i{r,+u^-rj-Uj)x{u,)xiuj). 

In Fourier space, the convolution products take a simple 
form so that 



5a;3(k) 



37rfi,r k'^5ap ^ kakj3 _ 
fcg k^-k^ + iO+ 



^X^(k) (15) 



where 5Q^(k) = J d^r e^^^'^ga/iir) is the Fourier trans- 
form of gai3 and x(k) is the one of x(i")- Second, the 
spontaneous emission rate F in Eqs. (|1I4I) has to be re- 
placed by 



r = rx'(ko 



(16) 



where ko is a vector of modulus equal to fco and of ar- 
bitrary direction. If one would treat the atomic motion 
quantum mechanically, as in [7[, for atoms trapped at 
the nodes of an optical lattice, x(u) — 0^(u) would 
be the probability distribution of the fluctuations u of 
the atomic position around a node r^, where </> is the 
underlying atomic center-of-mass wavefunction. Then 
Eq. (J14p would have a straightforward physical interpre- 
tation. Also f would simply be the elastic spontaneous 
emission rate, where the atomic center-of-mass after de- 
cay to the electronic ground state remained in the wave- 
function 0. In practice, a Gaussian choice for x is conve- 
nient, which corresponds to 



X(k) = e" 



(17) 



It is useful to know to which extent the results from the 
spatially smoothed model differ from the original model. 
For the Gaussian smoothing function, one then has the 
remarkable result that, when the width b is much smaller 
than all interatomic distances | r^ — r j | , one has the ap- 
proximate relation 



gani^i 



e ''"'''ffa/3(r»- 



(18) 



with an exponentially small error in 1/6^ [7[, that is one 
has the same Gaussian factor as for F. For the eigenvalue 
problem ([T]), this shows that the eigenvalues w — 17 of 
the spatially smoothed model may be related to the ones 
w — 47 of the original model by 



W — Wq — 17 : 



-fcnb 



(oj — Wo ~ il) 



(19) 



within an exponentially small error in 1/6^. For the 
steady state problem dU, it is found that the forced 
dipoles of the spatially smoothed model will (within an 
exponentially small error) coincide with the ones of the 
original model if one takes in the smoothed model the 
modified source frequency such that 



Wo 



= P. '^o" 



{uJs -Wo). 



(20) 



A. Density of states for the infinite periodic 
system 

In this subsection, we show how to recover Fourier 
space results of [^ for the density of states in the infi- 
nite periodic system, starting from the smoothed version 
of the real space Eq. ([T]). 

According to Bloch theorem, solutions of ([1]) can be 

taken of the form di_a ~ Oa e*'''^, where q is the Bloch 
vector, R is a vector of the Bravais lattice, the index 
/i labels primitive cells (for the diamond lattice, given 
by the combination of two shifted fee Bravais lattices, 
/i assumes two values), so that all atomic positions can 
be written as r^ = R -1- r'^^ where r^^^ is the position 
with respect to the Bravais lattice vector R. Injecting 
this ansatz in Eq. ([T|) modified according to Eqs. (|14ll6p . 
gives the eigenvalue problem 






(q)di' = hiuj - Wo - n) 



tM) 



(21) 



with 



■ afij^i/ 



(q) 



_ , , ^F ' 



5av 



+ 5]3„/3(R + r(^)-r(''))e-*i-^. (22) 
ReL 

Here indices a, /3 and /i, v label the direction and primi- 
tive cell, respectively, and eigenvalues w— 17 and eigenvec- 
tors (Ta depend on the choice of the cut-off smooth func- 
tion x(u), hence for the Gaussian choice as in Eq. (fT7| . 
they depend on the value of b. By considering the first 
contribution of Eq. iP^ . inside the square brackets, it is 
found from the inverse Fourier transform of (|15p that the 
tensor 30/3(0) is scalar (it is proportional to (5q^); fur- 
ther, using 1/{X + iO+) = P^ ~ *7r(5(X) and dS]), one 
finds that the imaginary part of gai3{0) exactly cancels 
with the F term. The second contribution, that is the 
sum over the Bravais lattice in (j22p . can be transformed 
with the Poisson summation formula. For the Gaussian 
smoothing function (IT71) . the real part of 5ct^(0) can be 
calculated explicitly; one obtains as in Q: 



^aii,l3u{Oi) = —Sa^S^i, 



l + 2(fco6)^ 

27ri/2(fcofe)3 



-Erfi(fco6)e 



-k^b^ 



-E 

Vt ^^ 



g,(K+q).(r(->-r<'^))|^^(K 



(23) 



where the wavevectors K run over the reciprocal lattice 
of the Bravais lattice, and Erfi is the imaginary error 
function. As expected for an infinite system, the matrix 
P is hermitian, so that 7 = 0. 

Turning back to the original problem ([T|) , that is in the 
absence of any smoothing function, we conclude for the 
infinite periodic system that the spectrum is real (7 = 0) 
and that huj — hujo is any of the eigenvalues of the matrix 



(q) 



lim ] 



(q), 



(24) 



as in the pcrturbativc limit of Q [that is for the eigen- 
frequencies close to Wq when Wp/wp — )- 0, tOp being the 
plasma frequency] . The resulting density of states is 



'"^'^sXil^^' 



^ - Wq,n) 



(25) 



where the integral over q is taken in the unit cell 2? = 
{^i=iQi^ii^\ — Qi < h} of the reciprocal lattice of 
basis (e.i)i<i<3, the sum over n runs over the all the 
eigenvectors of P(q) and Wq.n is the corresponding cigcn- 
frequency. 

For the Gaussian smoothing function, the limit of the 
band structure for 6 — >■ is computed in practice from 
the relation 



,,/3^(q) ~ PQ^,/3,y(q)e 



-kib' 



(26) 



which holds within an exponentially small error in 
(dmin/^)^ ^ 1 where dmin is the minimal interatomic 
distance [3, Q . Note that this relation, obtained for the 
particular case of a periodic system, is consistent with 
the general result P^ . and implies that the eigenvectors 
of P(q) essentially coincide with the ones of P(q). For the 
diamond, dmin = a-\/3/4. We used typically h = 0.05a, to 
which we applied the above & — ^ extrapolation formula 
to obtain the histogram in Figs. I1I3I5I 



B. Penetration depth for the infinite periodic 

system 

In this subsection we wish to derive, for an infinite 
system, the value of the penetration depth ^ and to con- 
firm that it depends on the considered direction of the 
direct space and that it diverges at the band edges, both 
properties having already been observed for a finite-size 
system in section UlI CI 

Hence, we have to solve Eq. (|4]) in presence of a forcing 



source dipole d^ 



placed in r^. The solutions 



we look for are the steady state dipole amplitudes di^a = 
'^R Q O'^ sach diamond lattice site of position r^ = R -|- 
Y^^\ where R belongs to the Bravais direct lattice. Since 
the scope is to determine the penetration length ^, we 
restrict ourselves to the case where the source frequency 
Ws is in the band gap. Then the dipole amplitudes are 
expected to decay exponentially at large distances, and 
one may introduce the Fourier transform 



<i - E 4".l^""-''- 



(27) 



ReL 



One applies this Fourier transform to the spatially 
smoothed version of Eq. @; for a Gaussian smoothing 
function, the source frequency is actually chosen to be uis 
given by Eq. (PO]) . which ensures that the forced dipole 
amplitudes arc essentially unaffected by the smoothing. 
In what follows, we can thus omit the bar (indicating the 



spatial smoothing) over the dipoles and the penetration 
depth. After calculations that closely resembles the ones 
of section IIV Al 



-E 






Vl 



^^„^(K + q)d). (28) 



KeRL 



One writes the formal solution of this linear system in 
terms of the inverse of the matrix P(q) — fi.(a)s — wg)!, 
where 1 is the identity; this inverse exists for all q since 
Ws is in the band gap of the spatially smoothed model. 
Then applying the inverse Fourier transform 



^R,Q 






(29) 



and using K • R = modulo 27r, one obtains the forced 
dipole amplitude on each lattice site: 



"r,q 



E E 

,3,7,iyKGRL 



(27r)3 



^gi(K+q)-(R+r('''-r, 



f[P(q)-ft(cS,-c.o)l] 'I ^;3^(K + q)J;. (30) 

A first application of Eq. ([50]) is to evaluate the dipole 
amplitudes from a numerical integration over q and, fit- 
ting them in a region of large values of R in some direc- 
tion u, to extract the penetration depth in that direction. 
Using up to 256'^ points in the numerical integration over 
q, this leads to the stars in Fig. [3 that compare well to 
the penetration depth extracted from the simulations on 
a finite size system in section fill CI Furthermore this ap- 
proach is numerically more efficient close to the borders of 
the band gap, where the penetration depth diverges and 
the finite size effects of the simulations become stronger. 

A second strategy to obtain the penetration depth from 
Eq. (pO| is to use the residue theorem. Since K + q spans 
R'^ when K spans the reciprocal lattice and q spans its 
unit cch V, and since P(q) = P(q + K), Eq. ^ can be 
rewritten as 



-'R.q 



E 



(2^ 



{[f{k)-h{Q,-ujo)t] U gpyik)df^. (31) 

To take the large R limit in the direction u, we set 

R = ru + 0(1) with r > 0. (32) 

We split the integration over k into an integral over the 
component fcii of k along u and over the transverse com- 
ponents k^ of k. Then k^ • R remains bounded, whereas 
u • R is divergent. 

First, we consider the integral over k\\ G M for a fixed 
kj^. The integrand involves the exponential factor e''^"''; 



since r > we close the integration contour with a half- 
circle (of diverging radius) in the upper complex plane 
[35| . Whereas the equation for A:||: 



<^fe||U+ki,?i — ^s 



(33) 



where a)k,n is the dispersion relation of the n-th band of 
eigenfrequencies for the spatially smoothed periodic sys- 
tem, has for sure no real solution since COg is in the band 
gap, it may have complex solutions k). with a positive 
imaginary part. Due to the occurrence of the inverse 
matrix involving P(k) in the integrand, such complex so- 
lutions provide poles in the half upper plane, which ac- 
cording to the residue theorem lead to the damped expo- 
nential exp(ifc|| r). If p3p admits several roots, or roots 
for various band index n, one has to keep the value uq 
of n and the root k.. leading to the smallest imaginary 
part, that provides the leading contribution in the large 
r limit. 

Then one has to remember that there is still an integral 
over kx, and that fc|| depends on k^. We thus face an 



integral of the form 






(Pk 



± a- "(ki)r 



(2^)2 



re II 



/(kj 



^fc|l'^fc[">(k^)u-hk^,„o 



(34) 



where the derivative of the band dispersion relation in 
the denominator originates from the residue of the pole 
in ku (k^) and the r-independent function / in the nu- 
merator is easily reconstructed from Eq. ((3T|) . To ob- 
tain an asymptotic equivalent of the integral ([34]) in the 
large-r limit, we use the saddle-point method: Eq. ([34]) 
is dominated by the contribution of the vicinity of the 



stationary point of the "phase", that is k 



(0) 



9k.fcr(kf) 



0. 



such that 



(35) 



,(0) 



As we shall see, in general k^"^ has complex coordinates 
(in the plane orthogonal to u) and one has to deform 
the integration domain of (IM]) to let the integration go 
through the stationary point [37| . Then one quadratizes 
the variation of the pole around the stationary point: 



u(0)/i,(0) 



k\{>{\i'l>+5\,. 



^^(kf) 



-5k_L-S5ki+0(5fci), (36) 



where the relevant deviations of k^ from the stationary 
point scale as l/r"^/"^. One finally gets the equivalent 



■,(0) 



i(M) 



e II 



')r 



/(kf) 



>oc afe„W^,,0)(j^(0))^^j^(0)^„^ 



(P5k 



i2n) 



■^ irSk±-BSkn 



(37) 
where the Gaussian integral provides a factor 1/r. The 
inverse of the penetration depth in direction u is thus 



k(u) 



1 



e(u) 



Im 



fcj^kf) 



(38) 



In general, this procedure is however difficult to use, 
even numerically, as one has to look for poles of the dis- 
persion relation for a wavevector k^*^) with three complex 
coordinates. An important and manageable limiting case 
is for a source frequency Ws very close to the lower bor- 
der Winf or the upper border uJsup of the band gap. The 
penetration depth is then expected to diverge, so that 
the imaginary components of the wavevector are small 
and its real components are close to the location qo in 
the Bloch vector space of the band gap border (such that 
^qo,„(, is equal to Wmf or (Dsup)- One can then quadratize 
the dispersion relation around the location of the border: 



Wqo+i5q,no 



dq-ASq+0{Sq^) 



(39) 



where A (rcsp. —A) is a positive definite matrix for the 
upper (resp. lower) border of the band gap. Note that, 
according to Eq. (|26p . A is related to its zero-6 limit A, 
that is to the matrix A of the original model, by 



-k:^b 



o" A 



(40) 



within an exponentially small error in 1/6^. Then the 
solution of ([33| obeying the stationarity condition ([35]) 
can be obtained analytically: 

k(°) ^ iti")(kf )u + kf ^ qo -I- ^.{u)^^ (41) 
with the expression for the inverse penetration depth 



k{u) = [(Wqo - UJs)u ■ A ^u] 



1/2 



(42) 



where A^^ is the inverse of the matrix A. In practice, 
one may find that the band gap border is obtained for 
several values of qo, due to symmetry properties (as it 
shall be the case for the diamond lattice) . At fixed direc- 
tion u, one then has to select the value of qo leading to 
the minimal value of «:(u) in Eq. (|1^ . Eqs. (J41I42|) are 
derived in the Appendix |X1 where the complete resulting 
expression for d^\^ is also given. 

A simple consequence of (|^^ is the asymptotic ex- 
pression for the maximal penetration depth at a given 
frequency w^, i-c. maximised over the direction u, close 
to a band gap border: 



in 



Ar. 



"boi-d \, CJbord ~ ^s 



1/2 



(43) 



where Amax is the eigenvalue of the matrix A of maximal 
modulus. 

We have explicitly evaluated the prediction (|42|) in the 
vicinity of the upper border of the band gap. Irrespective 
of the value of k^a, we find that the frequency Wsup of 
this upper border is reached on the so-called L point of 
the first Brillouin zone of the lattice, corresponding to 
qo = (ei +62-1- §3)72 = {■K/a){ej: + By + e^) [see Eq. (0 
for the values of the e^] , as it was already suspected in [8| . 
This point is so symmetric that all the six components of 
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the corresponding eigenvector of the matrix P are equal, 
which leads to the quite explicit expression 



r 



l^sup ^ I^Q — 



27rr 



l + 2(fco6)^ 
27ri/2(fco6)3 






Erfi {kab)e-''°''^ 



k^cV 



M 



)]- 



A-2 

ft.g 



j^/2 



(44) 



where K' = K + Qq. However, this frequency is also 
exactly reached for 13 other values of qo, so that 

Wqo = Wsup for 2qo G {±ei ± 63 ± eg, ±§1, ±63, ±63}. 

(45) 
For a given direction u, one thus calculates the 14 corre- 
sponding matrices A, which are all similar, and one keeps 
the one giving the smallest contribution to Eq. (|42|) . For 
u = ei and u = e^ this leads to the dashed line in 
the right part of Fig. [7^ and Fig. [71d respectively, in ex- 
cellent agreement with the numerical evaluation of (|30p 
and in good agreement with the finite-size simulations. 
Furthermore, for fcpa = 2 as in the simulations, the di- 
rection ei corresponds to the twice degenerate, maximal 
modulus eigenvalue A^nax of some of the 14 matrices A 
(the ones associated to qo = ±-^62 and qo = ±5^3) so 
that the maximal penetration depth ^max is obtained in 
that direction ei. Remarkably, for koa large enough (but 
smaller than the value kga ~ 5.14 leading to a closure of 
the gap), we find that the conclusion changes, and that 
the maximal penetration depth is now obtained in the 
direction (e^, -I- e^ -I- ez)/-\/3. This change suggests that 
there exists a magic value of koa such that the matrix A 
is scalar and, close to the upper bord of the band gap, 
the penetration depth is isotropic, which is confirmed by 
the diagonalisation of A that leads to [S^ : 



(koa) 



sup 
iso 



2.8632. 



(46) 



We have also explicitly evaluated the prediction of 
Eq. ((42|) in the vicinity of the lower border of the band 
gap. We have found that the frequency Winf of this lower 
border is obtained in 12 values qo of the Bloch vector, 
that weakly depend on k^a and that can be parameter- 
ized in terms of a single positive dimensionless unknown 
quantity a: 



Wqo =Winf for qoe{±cr(ei-e2),±cr(ei-e3), ±(7(62-63), 
± [a{ei + 62) + (2(7 - 1)63], ±[CT(ei + 63) + {2a - 1)62], 
±[a(e2+e3) + (2a-l)ei]} (47) 

where the basis vectors of the reciprocal of the fee lat- 
tice are given by Eq. ([6|) . Note that the last six elements 
of (|47p have a cr-independent component ±27r/a in the 
Cartesian basis, along e^, e.y, e^, respectively, and their 
components along the other two Cartesian axes are equal; 
these six elements are thus located on the straight line 
XU, where X and U are standard remarkable points of 
the first Brillouin zone of the diamond lattice. For the 



value koa = 2 taken in the figures, we numerically ob- 
tained a ~ 0.330 346. For those 12 values of qo, we have 
determined the 12 similar matrices A describing the lo- 
cal quadratization of cjq and we have kept, for a given u 
equal to ei or e^;, the one giving the smallest contribution 
to Eq. ([^ . This has led to the dashed line in the left 
part of Fig. [7^ and Fig. [Tjs respectively, again in excel- 
lent agreement with the numerical evaluation of (|30p and 
in good agreement with the finite-size simulations. For 
koa = 2, it is also found that ei is the eigenvector of two 
of the 12 similar A matrices [the ones corresponding to 
the last two elements of (|T7)) ] with the non degenerate, 
largest modulus eigenvalue ^max, so that the maximal 
penetration depth ^max is actually achieved in that di- 
rection, close to the lower border of the band gap. For 
larger values of fcoa, the situation can change to a maxi- 
mal penetration depth obtained along direction e^- This 
change occurs for the magic value 



inf 



ikoa)X 



2.9412 



(48) 



where a ~ 

value Amax 



0.353 740 and the maximal modulus eigen- 
of the matrices A is twice degenerate. 



C. States in the gap due to vacancies 

We now create a single vacancy in the periodic system 
(still using the spatially smoothed version), by removing 
the atom at the location Tig = Ro -I- r^^"^ that is at the 
lattice site Rq on the sublattice ^o- The eigenspectrum 
of the spatially smoothed version of ^ is expected to re- 
main real (7 = 0) but there may now be eigenvalues with 
id in the band gap of the periodic system, corresponding 
to states exponentially localized around the vacancy. As 
we will see, the corresponding w are given by Eq. (j53p . 

To look for such in-gap states, we use the following 
trick: Starting from a periodic system in presence of a 
source dipole in Ts (of imposed frequency w^ and ampli- 
tudes d^), we imagine that the vacancy on site r^g results 
from the coalescence of the corresponding forced dipole 
dio,a with the source dipole in the limit where the source 
location tends to the location of the vacancy: 



lim di, 



-di, Va. 



(49) 



In this case, the total dipole carried by the vacancy site 
vanishes, as if there was indeed a vacancy there. Ob- 
viously, condition P5|) can be satisfied only for specific 
values of w^ in the band gap of the spatially smoothed 
model, that we now determine. 

Writing Eq. ^ for R = Ro, /i = /lo, r^ 



Ro 



Jpo) 



j(A'o) 



and replacing d^"'^ with — d^, we obtain the homoge- 
neous linear system 



xQ^.,^;.o(q)d; (50) 



Pn,'^' 
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where we used (27r)'^ — VrlVl and wc called QQ^,/3,y(q) 
the non-scalar contribution to PQ^./3,y(q), that is the sec- 
ond contribution in the right-hand side of Eq. ([23]). In 
terms of matrices, 



P(q)=Al + Q(q), 



(51) 



where A is the coefhcient of the scalar contribution, that 
is of the first term in Eq. ^, A = -[^^^(O) + ihT/2] 
(this is independent of the direction a). We recognize 
a matrix product in Eq. (j50p . related to the sum over i/ 
and /3; wc then use 

[P(q)-Ml]-iQ(q) = t + {M - A)[¥{q) - hSt]-\ (52) 



with 5 — (Ds — wq- The contribution to (|50p of 1 in 
that expression exactly reproduces the term d^ of the 
left-hand side of ((50| . since the integral over q on the 
primitive cell V of the reciprocal lattice is equal to Vrl- 
Simplifying the remaining contribution by the factor M— 
A. it remains 

= / ^ y {[P(q) - K^s - c^o)i]-'} <, Va. 

J-D ^RL U v^^ V "^ J JQ/io,7^Q 7' 

(53) 
This must have a non-zero solution for the source dipole, 
which is equivalent to requiring that the <Ds-dependent 3x 
3 hcrmitian matrix in (j53|) has a zero eigenvalue. To show 
that the condition ([55]) is not only sufficient, but also 
necessary, we have performed an alternative calculation, 
presented in Appendix [B] that has also the advantage of 
including the case of several vacancies. 

For the diamond lattice, we have evaluated numerically 
the integral over the Bloch vector q in Eq. (|53p . We 
then find that the resulting 3x3 hcrmitian matrix is 
scalar. As the eigenvalues of that matrix are increasing 
functions of ujg, as can be shown with the Hellmann- 
Fcynman theorem, this implies that there is at most one 
solution for uj^ in the band gap. Numerically, we find that 
there is a solution, whose value [after extrapolation to 
6 — >■ using Eq. p9)) ] for fego = 2 is indicated by a vertical 
dotted line in Fig. [Sj in agreement with a peak location 
in the density of states in the numerical simulations. 

In the case of several vacancies, we can extend our anal- 
ysis as described in appendix [BJ By numerical solution of 
Eq. (jB7p . we have then investigated the in-gap states for 
two vacancies on sites separated by R2 — Ri = 0, ei or 
ae^;, being either on the same sublattice (/ii = /i.2) or on 
different sublattices {jli ^ /i2)- In most cases, we have 
found allowed frequencies close to the one of the single- 
vacancy state, within the width of the central peak in the 
inset of Fig. [31 those states can not be resolved in that 
figure and we have not indicated them. For the two ge- 
ometries specified in the caption of Fig. [3l we have found 
frequencies of two-vacancy states that are clearly out of 
the central peak, see the red and blue vertical dotted 
lines; in particular, the prediction with (w — cjo)/r ~ —4 
seems to match quite well the very clear secondary peak 
that emerges in the figure for increasing concentration of 
vacancies. 



V. CONCLUSION 

Three-dimensional periodic arrangements of extended 
scattering objects leading to an omnidirectional band gap 
for light have been known since the 90's, starting from the 
diamond lattice configuration of dielectric microspheres 
of [30| . In the case of a periodic ensemble of point-like 
scatterers, the technical issues affecting the calculation of 
the band structure of light have been solved only recently 
[3, 1^1 13 1 which has allowed to show that the diamond 
lattice can also lead to a photonic band gap in the point- 
like case [i]. 

With cold atom experiments, a diamond-like ensemble 
of point-like scatterers is in principle realizable, provided 
that oneproduces, in the appropriate optical lattice ge- 
ometry [3, UM ' ^ high quality Mott phase of atoms [5|, |^| 
having an optical transition between a spin zero ground 
state and a spin one electronic excited state |3l| . In prac- 
tical realizations, there will be of course unavoidable de- 
viations from the ideal infinite periodic case, that we have 
quantified in the present work with numerical solutions 
of linearly coupled dipoles equations with about 3 x 10''^ 
particles. 

A first issue is due to effects of the finite size of the 
atomic medium. Rather than having a band structure, 
light has a continuous spectrum of scattering states; by 
analytic continuation to the lower half of the complex 
plane, however, it is more physical to consider, as we have 
done, the discrete complex eigenfrequencies w — 17 of the 
resonances of the system. In the distribution function of 
a;, the forbidden gap remains visible in our simulations. 
It remains actually quite visible if one restricts to the 
resonances with a half decay rate 7 much smaller than 
the free space single atom spontaneous emission rate F; 
such a filtering of the resonances could be realized exper- 
imentally by performing a frequency measurement after 
an adjustable time delay, during which the short-lived 
resonances decay and are suppressed. Amusingly, a nar- 
row peak in the distribution function of w was observed 
close to the center of the infinite system band gap, when 
the finite size atomic medium has a cubic shape; such a 
peak, absent when the medium has a spherical shape, is 
a very clear finite size effect. 

A second issue is due to vacancies inside the atomic 
medium. For a concentration of a few per cent of vacan- 
cies, narrow peaks emerge in the distribution function 
of uj inside the gap. We were able to identify several 
of these peaks as corresponding to the frequencies of lo- 
calised states around one or two close vacancies in an 
otherwise infinite periodic medium. At higher concen- 
trations of vacancies, e.g. 20%, with no filtering on 7, 
the gap disappears. 

From our finite size sample, we have shown that one 
can quite accurately extract the penetration depth ^ of 
the light in the medium, and that the obtained values 
compare well with independent calculations in a periodic 
medium. Away from the borders of the band gap, ^ as 
a function of the imposed field frequency 0;^ exhibits a 
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plateau at a remarkably low value, between 0.5a and a, 
where a is the lattice constant of the underlying fee lat- 
tice. Close to the borders Wbord of the band gap, one can 
even directly observe, in our finite size system, the onset 
of the divergence of ^ as 1/\lOs — Wbordl^^^, with a pref- 
actor close to our analytical predictions. We have also 
observed from the simulations that ^ is anisotropic (it de- 
pends on the direction of space), in agreement with our 
theoretical analysis, and that this anisotropy becomes 
quite pronounced close to the lower border of the band 
gap. 
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Appendix A: Penetration depth 

In this Appendix, for the spatially smoothed model, 
we derive the results (|41l42p for the penetration depth 
in the direction u at a frequency Cjs close to a border of 
the band gap, which justifies the use of the quadratizcd 
dispersion relation ([M|) around the Bloch vector qo, and 
we give the large-distance equivalent of the forced dipole 
amplitude, as obtained from the saddle-point method. 

As short-hand notations, we introduce z = fcy — 9o|l f-nd 
X = k^ — qo_L as the components along u and in the plane 
orthogonal to u of the vector k — Qq. We also introduce 
the frequency deviation from the nearest band border, 
A = ajq„ — LOs- Then Eq. ([55]) reduces to a degree- two 
equation for z: 



z^n ■ Au. + 2zx • v4u + X • ^x -I- A = 



(Al) 



Furthermore, z has to be stationary with respect to a 
variation of k^, see Eq. (|35p . Differentiating the trino- 
mial (jAip with respect to x, and using 9x2 = 0, one 
obtains the vectorial equation zQ Am + QAQyi ~ where 
Q projects orthogonally to u. The solution is 



x = ~z{QAQ)-^A\i 



(A2) 



where the matrix inverse is intended within 
the vectorial plane orthogonal to u. Insert- 

ing this solution into Eq. (|A1|) and using 
[PAP - PAQ[QAQ)-^QAP]PA-^P = P where 
P = 1 — Q is the orthogonal projector on u (see relation 
(B.23) of §III.B.2 in ^), one obtains 



z = zk(u) with k(u) given by Eq. 



(A3) 



Similarly, injecting the closure relation P + Q = 
1, _ one _finds ^[u - {QAQ)-^Au] = [PAP - 
PAQ {QAQ)~^ QAP]u 
This gives as in Eq. (HI 



(PA-ip)-iu = u/(u • A-^u). 



u-{QAQ)-^Au 



u 



A- 



(A4) 



To determine the residue appearing in (j37|) . one takes 
the derivative of the trinomial (|Aip with respect to z for 
a fixed X. Using the previous relations one obtains 



d, 



k«^Wm,no 



2zu ■ Au + 2yL- Ax = 



2iK(u) 
(u-^-iu)' 



(A5) 



Next, we determine the matrix B in Eq. p7p originating 
from the quadratization of z around the stationary point 
x. A first order variation Sx induces a second order vari- 
ation 6z. Performing these variations in Eq. (jAip up to 
second order in Sx and up to first order in Sz, and using 
the previous relations, we obtain 



B 



2k(u) 



(u- A-^u)QAQ. 



(A6) 



We conclude that the matrix iB appearing in the Gaus- 
sian integral ([STjl is negative, which justifies the fact that 
the saddle point is approached along the real axis direc- 
tion as in (j37p . If one performs the Gaussian integral, 
Eq. ([37|) reduces to 



4"^ 

tt.a 



-K.(u)r 



/(kf) 



Ainr 



[dct(QAQ)]-i/2. (A7) 



The determinant in that expression is conveniently trans- 
formed as det{QAQ) = (u-^~^u) det A using the expres- 
sion of the matrix of A~^ in terms of the comatrix of A 
(in an orthonormal basis containing the direction u). 

To obtain our final asymptotic form for the forced 
dipole amplitude, we note that, for any acceptable vec- 
tor k^"^ of the pole plus saddle-point analysis, k^°^ -I- K 
is again acceptable, where K is any vector of the recipro- 
cal lattice; this is due to the periodicity of the dispersion 
relation Wk.no- We also include a sum over possibly de- 
generate Bloch vector qg leading to the same value uiqg^no 
(as discussed in the main text). We also note that, when 
R -^ 4-0O, 

||(i-iA)i/2(R+r('')-r,)|| = -z(k(o)-qo)-(R+r('^)-r,) 

+ 0(1) (A8) 

which gives a simple physical interpretation to the ex- 
pression (PT|) of k(°^ : The apparently obscure correction 
to qo in (|4ip simply originates from the fact that what 
more precisely matters in the asymptotic behavior of the 
dipole amplitudes is not ru but really the vectorial dis- 
tance R-f-r'"^' — Ts between the considered lattice site and 
the source. Finally, we obtain, for Wg close to a border 
of the band gap, the asymptotic equivalent for R — >■ cxd: 






E 

qo 



/ A \ i/2g»qo-R.e-ll(^"'^)'^'(^+'-''' 



\detAj 47r;i||(y4-iA)i/2(R + r('^) -rs 



X E E e'^'>°+^)-e^'"'--'e;V?;^*^^.(qo+K)d;, 

/J./j.i/KeRL 

(A9) 
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where ^L'/j" are the components of the normaUzed eigen- 
vector of P(qo) of eigenvalue ojqo,noj we approximated 
k'°) by qo in the argument of gp^, and the square root 
[A~^Ay/'^ of the matrix A~^A is well defined since this 
matrix is positive. Note that the second line of (|A9|) does 
not depend on R. 



where S is the Kronccker symbol and we have introduced 
the auxiliary unknowns 



^ g(R, + r^. - R' - v^,)d'^?. (B4) 



Appendix B: A general vacancy calculation 



We consider here the infinite periodic system, with a 
finite number of vacancies at nodes (Ri,^^), 1 < i < n, 
where we recall that R belongs to the fee Bravais lattice 
and /x labels the sublattices. The scope is to determine 
the frequencies to of the localised states that can exist, 
due to the presence of the vacancies, in the band gap of 
the periodic system, in the spatially smoothed version of 
the model. 

The idea is to formally introduce, in the coupled equa- 
tions for the dipoles, fictitious dipoles carried by the 
vacancies. Among the physical dipoles, the spatially 
smoothed version of Eq. ([T]) holds: 



Then taking the Fourier transform (P7)) of Eq. (jB3p and 
using 



p,' i—l 

Since the frequency w is in the gap, the matrix is in- 
vertible, and taking the inverse Fourier transform, one 
obtains 



= (A - M)d*^' 






-r^-R'-V)<' 



(p') 



(Bl) 



Here the prime over the summation symbol means that 
the sum is restricted to the physical dipoles, d = ui — luq is 
the detuning from the atomic resonance, A is defined be- 
low Eq. (j5ip , and we have used for conciseness an implicit 
vectorial notation for the dipoles and an implicit matrix 
notation for g. For the fictitious dipoles, the equation is 
that they are equal to zero: 



= (A 



hS)df^ 



Vie{l,...,n}. 



(B2) 



This allows to formally extend the sum in Eq. (jBTj) to 
the fictitious dipoles, that is one can remove the prime 
over the summation symbol. One can then merge the 
two series of equations using the usual plus-minus trick: 
for all R in the Bravais lattice and for all sublattices /i, 
one requires that 

= (A - M)4^ + ^ 5(R + i-p - R' - r,.')4' ^ 



R-'.m' 



XI Vr/^^/^'*' (^^) 



i=l 



n p 



i^piq.(R-R.; 



Vrl 



{[P(q) - Ml] -1 },,,,, s. 



(B6) 



Expressing the fact that the fictitious dipoles are all 
equal to zero, we find the homogeneous system of equa- 
tions: 



E 



■D 






^.q.(R,-R.){[p(q) _ ^^jj-l| ^^-^ ^ 0^ 



(B7) 

to be satisfied Vj S {l,...,ri}. The acceptable in- 
gap frequencies are such that the system admits a non- 
identically zero solution (si)i<i<„, that is the determi- 
nant of the corresponding in x 3n matrix must vanish. 
In the case of a single vacancy, this reproduces Eq. ([55]) . 



Finally we have performed the consistency check that, 
if one replaces in Eq. (|B4p the dipoles in terms of the 
auxiliary unknowns Sj, as given by (jB6p . one recovers 
exactly the same system as (|B7p . using Eqs. (|22l5ip and 
the fact that the integral over q on the primitive cell T) 
of the reciprocal lattice is equal to its volume Vrl ■ 



[1] G. Grosso, G. Pastori-Parravicini, Solid State Physics 

(Academic Press, 2000). 
[2] A.M. Afanas'ev and Yu. Kagan, Sov. Phys. JETP 25, 

124 (1967); G.B. Smirnov, Y.V. Shvydko, JETP Letters 

35, 505 (1982). 



[3] J.J. Hopfield, Phys. Rev. 112, 1555 (1958); V. Agra- 
novich, Sov. Phys. JETP 37, 307 (1960). 

[4] J.D. Joannopoulos, S.G. Johnson, J.N. Winn, and R.D. 
Meade, Photonic Crystals: Molding the Flow of Light 
[Princeton University Press, Princeton, NJ, 2008] (2nd 



14 



Edition); F. Zolla, G. Renversez, A. Nicolet, B. Kuhlmey, [25] 

S. Guenneau, D. Felbacq, A. Argyros, and S. Leon-Saval, [26] 

Foundations of Photonic Crystal Fibres [Imperial College 

Press, London, 2012] (2nd Edition). 
[5] M. Greiner, O. Mandel, T. Esslinger, T.W. Hansch, and 

I. Bloch, Nature 415, 39 (2002). 
[6] M. Anderlini, P.J. Lee, B.L. Brown, J. Sebby-Strabley, 

W.D. Phillips, J.V. Porto, Nature 448, 452 (2007). 
[7] M. Antezza and Y. Castin, Phys. Rev. Lett. 103, 123903 

(2009). [27 

[8] M. Antezza and Y. Castin, Phys. Rev. A 80, 013816 

(2009). 
[9] Masao Takamoto, Feng-Lei Hong, Ryoichi Higashi, [28 

Hidetoshi Katori, Nature 435, 321 (2005). 
[10] T.L. Nicholson, M.J. Martin, J.R. Wilhams, B.J. Bloom, [29 

M. Bishof, M.D. Swallows, S.L. CampbeU, and J. Ye, 

Phys. Rev. Lett. 109, 230801 (2012) , and references 

therein. [30 

[11] A.D. Ludlow, T. Zelevinsky, G.K. Campbell, S. Blatt, 

M.M. Boyd, M.H.G. de Miranda, M.J. Martin, J.W. [31 

Thomsen, S.M. Foreman, Jun Ye, T.M. Fortier, J.E. Stal- 

naker, S.A. Diddams, Y. Le Coq, Z. W. Barber, N. Poh, [32 

N.D. Lemke, K.M. Beck, and C.W. Gates, Science 139, 

1805 (2008). 
[12] D.V. van Coevorden, R. Sprik, A. Tip, and A. Lagendijk, 

Phys. Rev. Lett. 77, 2412 (1996). [33 

[13] P. de Vries, D.V. van Coevorden, A. Lagendijk, Rev. 

Mod. Phys. 70, 447 (1998). [34 

[14] J. A. Klugkist, M. Mostovoy, and J. Knoester, Phys. Rev. 

Lett. 96, 163903 (2006). 
[15] A. Schilke, C. Zimmermann, P. W. Courteille, and W. 

Guerin Phys. Rev. Lett. 106, 223903 (2011); A. Schilke, 

C. Zimmermann, and W. Guerin, Phys. Rev. A 86, [35 

023809 (2012). 
[16] A. Schilke, C. Zimmermann, P. W. Courteille, and W. 

Guerin, Nature Photonics 6, 101 (2011) [36 

[17] S. Rist, C. Menotti, and G. Morigi, Phys. Rev. A 81, 

013404 (2010). rg^ 

[18] I. Carusotto, M. Antezza, F. Bariani, S. De Liberato, and 

C. Giuti Phys. Rev. A 77, 063621 (2008). 
[19] H. Zoubi, H. Ritsch, Phys. Rev. A 76, 013817 (2007). 
[20] O. Morice, Y. Castin and J. Dahbard, Phys. Rev. A 51, [38 

3896 (1995). 
[21] D. Felbacq, and M. Antezza, SPIE Newsroom (2012) 

[DOI: 10.1117/2.1201206.004296], and references therein. 
[22] U. Fano, Phys. Rev. 103, 1202 (1956). 
[23] A. Chelnokov, S. Rowson, J.-M. Lourtioz, V. Berger, J.- 

Y. Courtois, J. Opt. A: Pure Appl. Opt. 1, L3 (1999). 
[24] O. Toader, T.Y. Chan, and S. John, Phys. Rev. Lett. 92, 

043905 (2004). 



D. Yu, Phys. Rev. A 84, 043833 (2011). 
It is worth stressing that the often used scalar model for 
light has the evident advantage of drastically reducing the 
numerical effort, but also the disadvantage of providing 
a qualitatively and quantitatively wrong description of 
the physical system. For instance, it is possible to show 
that already for a simple cubic atomic lattice, the scalar 
model, in contradiction with the vectorial one, predicts 
the presence of a band gap. 

Y. Bidel, B. Klappauf, J.C. Bernard, D. Delande, G. 
Labeyrie, C. Miniatura, D. Wilkowski and R. Kaiser, 
Phys. Rev. Lett. 88, 203902 (2002). 

J.D. Jackson, Classical Electrodynamics, 2nd ed. (Wiley, 
New York, 1975). 

C. Cohen- Tannoudji, J. Dupont-Roc, G. Grynberg, Pro- 
cessus d 'interaction entre photons et atomes, InterEdi- 
tions/Editions du CNRS (Paris, 1988). 
K.M. Ho, C.T. Chan, CM. Soukouhs, Phys. Rev. Lett. 
65, 3152 (1990). 

T. Fukuhara, S. Sugawa, M. Sugimoto, S. Taie, Y. Taka- 
hashi, Phys. Rev. A 79, 041604 (2009). 
Although an optical diamond lattice was to our knowl- 
edge not realized yet in the lab, the technique to be ap- 
plied, elaborating on the ideas of [23] , is perfectly known 

mm 

Our convention amounts to omitting the S{r) term in 
Eq. (3) of [13]. ^^ 

Note that, according to Eqs. (|19I20[) to come, this ra- 
tional fraction of the source frequency is the same for 
the original model and the Gaussian spatially smoothed 

model, (tJs — cDinf)/(cDsup— tJs) ~ {lOs — UJini)/{uJBup — UJs), 

within an exponentially small error in l/&^. 
To this end, the Gaussian smoothing function x is not 
appropriate. One can rather take x(r) oc e~^'''/r, whose 
Fourier transform is a Lorentzian. 

If there are several stationary points, one has to keep the 
one leading to the smallest imaginary part of fc:, . 
This is why the naive minimization of the imaginary part 
of fc|| (kx) over real-component kx only gives an upper 
bound on the penetration length in the direction u. 
For koa below that value, Amax is twice degenerate and 
the corresponding eigenspace is the plane orthogonal to 
(e^: + By + e^)/^ for qo = (ei + §2 + e3)/2 and the 
plane orthogonal to (— e^; + By + e,)/y/3 for qo — ei/2. 
For koa above that value, Amax is not degenerate and the 
corresponding eigenvector is (e^ -I- e^ + ez)/V3 for qo = 
(ei +e2+e3)/2, and {—ex + ey + ez)/VSior qo = ei/2. 



